PHASE RETRIEVAL FOR IMAGING PROBLEMS. 

FAJWEL FOGEL, IRENE WALDSPURGER, AND ALEXANDRE D' ASPREMONT 

Abstract. We study convex relaxation algorithms for phase retrieval on imaging problems. We show that 
structural assumptions on the signal and the observations, such as sparsity, smoothness or positivity, can be 
exploited to both speed-up convergence and improve recovery performance. We detail experimental results in 
^^ molecular imaging problems simulated from PDB data. 

o 

(N 

p^ 1. Introduction 

"^ Phase retrieval seeks to reconstruct a complex signal, given a number of observations on the magnitude 

0^ of linear measurements, i.e. solve 

^ find X 



U 
O 



such that \Ax\ 



in the variable x G C^, where A ^ MJ^^p and 6 G M^. This problem has direct applications in X-ray 

and crystallography imaging, diffraction imaging, Fourier optics or microscopy for example, in problems 

4^ where physical limitations mean detectors usually capture the intensity of observations but cannot recover 

C^ their phase. In what follows, we will focus on problems arising in diffraction imaging, where A is usually 

2 a Fourier transform, often composed with one or multiple masks (a technique sometimes called ptychogra- 

■"^ phy). The Fourier structure, through the FFT, often considerably speeds up basic linear operations, which 

^_^ allows us to solve large scale convex relaxations on realistically large imaging problems. We will also ob- 

^ serve that in most of the imaging problems we consider, the Fourier transform is very sparse, with known 

^ support (we lose the phase but observe the magnitude of Fourier coefficients), which allows us to consider- 

j.^ ably reduce the size of our convex phase retrieval relaxations. 

f"^ Because the phase constraint \Ax\ = 6 is nonconvex, the phase recovery problem (1) is non-convex. 

^ Several greedy algorithms have been developed (see [Gerchberg and Saxton, 1972; Fienup, 1982; Griffin 

O and Lim, 1984; Bauschke et al., 2002] among others), which alternate projections on the range of A and 

^^ on the nonconvex set of vectors y such that |y| = \Ax\. While empirical performance is often good, these 

algorithms can stall in local minima. A convex relaxation was introduced in [Chai et al., 201 1] and [Candes 

.^ et al., 2011b] (who call it PhaseLift) by observing that \Ax\'^ is a linear function of X = xx* which is a 

S/ rank one Hermitian matrix, using the classical lifting argument for nonconvex quadratic programs developed 

H in [Shor, 1987; Lovasz and Schrijver, 1991]. The recovery of x is thus expressed as a rank minimization 

problem over positive semidefinite Hermitian matrices X satisfying some linear conditions, i.e. a matrix 

completion problem. This last problem is approximated by a semidefinite program which has been shown to 

recover x for several (random) classes of operators A [Candes et al., 2011b,a]. More recently, [Waldspurger 

et al., 2012] showed that the phase retrieval problem (1) can be written as an extension of the MAXCUT 

combinatorial graph partitioning problem over the unit complex torus, allowing fast algorithm for solving 

semidefinite relaxations of MAXCUT to be applied to phase retrieval. 

On the experimental side, phase recovery is classical problem in Fourier optics for example [Goodman, 
2008], where a diffraction medium takes the place of a lens. This has a direct applications in X-ray and 
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crystallography imaging, diffraction imaging or microscopy [Harrison, 1993; Bunk et al., 2007; Johnson 
et al., 2008; Miao et al, 2008; Dierolf et al, 2010]. 

Here, we implement and study convex relaxation algorithms for phase retrieval on imaging problem 
instances where A is based on a Fourier operator. We show how structural assumptions on the signal and 
the observations (e.g. sparsity, smoothness, positivity, support, etc.) can be exploited to both speed-up 
convergence and improve recovery performance. The paper is organized as follows. Section 2 briefly recalls 
the structure of some key algorithms used in phase retrieval. Section 3 describes applications to imaging 
problems and how structural assumptions can considerably reduce the cost of solving large-scale instances.. 
Section 4 details some numerical experiments while Section 5 describes the interface to the numerical library 
developed for these problems. 

Notations. We write Sp (resp. H^) the cone of symmetric (resp. Hermitian) matrices of dimension p ; S^ 
(resp. H^) denotes the set of positive symmetric (resp. Hermitian) matrices. We write || • ||p the Schattenp- 
norm of a matrix, that is the p-norm of the vector of its eigenvalues (in particular, 1 1 • 1 1 oo is the spectral norm). 
We write A^ the (Moore-Penrose) pseudoinverse of a matrix A, and AoB the Hadamard (or componentwise) 
product of the matrices A and B. For x G MP, we write diag(x) the matrix with diagonal x. When X G Hp 
however, diag(X) is the vector containing the diagonal elements of X. For X G H^, X* is the Hermitian 
transpose of X, with X* = (X)^. Finally, we write 6^ the vector with components bf,i = 1, . . . , n. 

2. Algorithms 

The phase problem (1) can be written explicitly in terms of the phase variable instead of the signal x. In 
the noiseless case, we can write the constraint \Ax\ = 6 as Ax = diag(6)ix, where i/ G C^ is a phase vector, 
satisfying l^x^l = 1 for i = 1, . . . , n, so problem (1) becomes 

minimize H^x — diag(6)i^||2 
subject to \ui\ = 1 

where we optimize over both phase u ^ C^ and signal x G C^. While the objective of this last problem is 
jointly convex in (x, u), the phase constraint l^x^l = 1 is not. This last problem can also be written 

minimize llAx — tvIIo 

I (3) 

subject to \y\ ^ b 

where we now optimize over both phased observations y G C^ and signal x G C^. Several greedy algorithms 
attempt to solve this last problem using variants of alternating projections and we detail some of the most 
classical examples in the paragraphs that follow. 

2.1. Greedy algorithms. Several greedy algorithms based on alternate projections have been developed to 
solve (1). The algorithm Gerchberg-Saxton by [Gerchberg and Saxton, 1972] for example seeks to recon- 
struct y = Ax and alternates between orthogonal projections on the range of A and normalization of the 
magnitudes |^| to match the observations b. The cost per iteration of this method is minimal but convergence 
(when it happens) is often slow. 

Algorithm 1 Gerchberg-Saxton. 



Input: An initial ; 


t/i G F, i.e. 


such that y = 


b. 


1: for fc = 1,... 


, AT - 1 do 






2: Set 


y^' ^ 


' liAA^y'hl 


, 


3: end for 








Output: yN e F. 









1, . . . , n. (Gerchberg-Saxton) 



A classical "input-output" variant, detailed here as algorithm Fienup, introduced by [Fienup, 1982] adds 
an extra penalization step which usually speeds up convergence and improves recovery performance when 
additional information is available on the support of the signal. Of course, in all these cases, convergence to 
a global optimum cannot be guaranteed but empirical performance is often quite good. 

Algorithm 2 Fienup 

Input: An initial y^ G F, i.e. such that |y^| = 6, a parameter /3 > 0. 
1: for A: = 1,...,7V- Ido 
2: Set 

_ {AA^y% ._, 



3: Set 

4: end for 
Output: yN eF. 



v\^^ = v\ - ^{v\ - h^i) (Fienup) 



2.2. Semidefinite programming. In [Chai et al., 2011; Candes et al., 2011b], the phase recovery problem 

is then written as 

minimize Rank(X) 

subject to Tr (aia*X) = 6f , z = 1, . . . , n 

X^O 

in the variable X G Hp, where X = xx* when exact recovery occurs. This last problem can be relaxed as 

minimize Tr(X) 

subject to Tr (aia*X) = 6^, z = 1, . . . , n (PhaseLift) 

X ^0 

which is a semidefinite program (called PhaseLift by Candes et al. [2011b]) in the variable X G Hp. This 
problem is solved in [Candes et al., 2011b] using first order methods. 

In [Waldspurger et al., 2012], it was shown that the phase problem retrieval problem (1) could be split 
between phase and signal reconstruction, exploiting the fact that given the phase, signal reconstruction is 
a simple least squares problem. Problem (2) is non-convex in the phase variable u. However, given u we 
obtain x as 

x^ J^ diag(6)i/ (4) 

where J^ is the pseudo inverse of A. Replacing x by its value in (2), the phase recovery problem becomes 

minimize u'Mu 

subject to \ui\ — \^ z = l,...n, 

in the variable i/ G C^, where the Hermitian matrix 

M = diag(6)(I - AA^) diag(6) 

is positive semidefinite. [Waldspurger et al., 2012] detailed greedy algorithm Greedy to optimize (5) in the 

phase variable. 

A convex relaxation to (1) was also derived in Waldspurger et al. [2012] using the classical lifting argument 

for nonconvex quadratic programs developed in [Shor, 1987; Lovasz and Schrijver, 1991]. This relaxation 

is written 

min. Tr(f/M) rphaseCurt 

subject to diag(f/) = 1, f/ ^ 0, ^i-nasecutj 

3 



Algorithm 3 Greedy algorithm in phase. 



Input: An initial u e C^ such that l^x^l = 1, z = 1, . . . , n. An integer A^ > 1. 



1: for A: = 1,.. 


,A^do 


2: for z = 1, . 


. .ndo 


3: Set 





Ui = ^'^' ' ' (Greedy) 



T.J^^MjiU 



4: end for 
5: end for 

Output: u ^£P such that l^/^l = 1, z = 1, . . . , n. 



which is a semidefinite program (SDP) in the matrix t/ G H^. This problem has a structure similar to the 
classical MAXCUT relaxation and instances of reasonable size can be solved using specialized implemen- 
tations of interior point methods. [Waldspurger et al., 2012] solve larger instances of this problem using the 
block-coordinate descent algorithm PhaseCut. 

Algorithm 4 Block Coordinate Descent Algorithm for PhaseCut. 

Input: An initial U^ = I^ and v > Q (typically small). An integer A^ > 1. 
1: for A: = l,...,A^do 

2: Picki E [l,n]. 
3: Compute 

u = Ufc icMic^i and 7 = u*Mic^i (PhaseCut) 



4: If 7 > 0, set 



else 



m'-ut^'*--/-^- 



Tj-k+l _ Tj-k+l^ 



5: end for 

Output: A matrix 1/^0 with diag(t/) = 1. 



3. Imaging PROBLEMS 

In the problems we study here, the matrix A is usually structured and we have significant structural 
information on both the signal we seek to reconstruct (regularity, etc.) and the observations (power law 
decay in frequency domain, etc.). Many of these additional structural hints can be used to both speedup 
convergence and improve phase retrieval performance. The paragraphs that follow explore these points in 
more detail. 

3.1. Exploiting Fourier structure. In practical applications, because of the structure of the linear opera- 
tor A, we may often save computational time by computing quickly the required matrix-vector products. 
We explain the case where A corresponds to the Fourier transform with masks. Let /i, ..., /^^^ G C^ be the 
illumination masks. The image by A of some signal x G C^ is 

/^(/lox)\ 



Ax 



[jXh o x)J 



and the pseudo-inverse of A also has a simple structure: 

(yi\ k 

where li is the dual filter of // which is 



1=1 



''-'■'[V'-n- 

With the Fast Fourier Transform, the application of A or A'^ only requires 0{kp\og(p)) operations. So for 
any v G C^, Mv = diag(6)(I — AA'^) di3.g{b)v may be calculated in 0{kplog{p)) operations instead of 
0{k^p^) for naive matrix-vector multiplications. 

In algorithms Greedy and PhaseCut, we also need to extract quickly columns from M without having 
to store the whole matrix. Extracting the column corresponding to index i in block / < k redues to the 
computation of AA^i^i where 5i^i G C^^ is the vector whose coordinates are all zero, except the z-th one of 
/-th block. If we write 5i G C^ the Dirac in z, the preceding formulas yields 

A^^(/lo//)\ 

AA^.^i = : 

\5i^:F{IkoI[)) 

Convolution with 5i is only a shift and vectors J^{Is o I'l) may be precomputed so this operation is very fast. 

3.2. Low rank iterations. In instances where exact recovery occurs, the solution to the semidefinite pro- 
gramming relaxation (PhaseCut) has rank one. It is also likely to have to have low rank in a neighborhood 
of the optimum. This means that we can store a compressed version of the iterates U in algorithm PhaseCut 
in the form of their SVD VV where V G C^^^. Each iteration updates a single row/column of U which 
corresponds to a rank two update of U, hence updating the SVD means computing a few leading eigenvalues 
of the matrix VV + LL"" where L G C^^^. This update can be performed using Lanczos type algorithms 
and has complexity 0{kn log n). Compressed storage of U saves memory and also speeds-up the evaluation 
of the vector matrix product Uic^icMic^i which costs 0{nk) given a decomposition Uic^ic = VV, instead of 
0{in?) using a generic representation of the matrix U. 

3.3. Bounded support. In many inverse problems the signal we are seeking to reconstruct is known to be 
sparse in some basis and exploiting this structural information explicitly usually improves signal recovery 
performance. This is for example the basis of compressed sensing where li penalties encourage sparsity 
and provide recovery guarantees when the true signal is actually sparse. 

The situation is a lot simpler in some of the molecular imaging problems we are studying below since the 
electron density we are trying to recover is smooth, which means that its Fourier transform will be sparse, 
with known support. While we lose the phase, we do observe the magnitude of the Fourier coefficients. This 
allows us to considerably reduce the size of the SDP relaxation without losing much reconstruction fidelity, 
i.e. in many cases we know that a significant fraction of the coefficients of h are close to zero. From a 
computational point of view, sparsity in h allows us to solve a truncated semidefinite relaxation (PhaseCut). 

Indeed, without loss of generality, we can reorder the observations such that h = (6f , 0)^. Similarly, we 
note 

Using the fact that 62 = 0, the matrix M in the objective of (5) can itself be written in blocks, that is 






Figure 1. Electronic density for the caffeine molecule (left), its 2D FFT transform 
(diffraction pattern, center), the density reconstructed using 2% of the coefficients with 
largest magnitude in the FFT (right). 



Since 62 = 0, any complex vector with coefficients of magnitude one can be taken for U2 and the optimiza- 
tion problem (5) is equivalent to 



(6) 



minimize ulMiUi 

subject to \ui. 1 = 1, i = 1, . . . n, 

in the variable ui G C^^ , where the Hermitian matrix 

Ml = diag(6i)(I - Ai(At)i) diag(6i) 

is positive semidefinite. This problem can in turn be relaxed into a PhaseCut problem which of course is 
usually considerably smaller than the original (PhaseCut) problem. 

3.4. Positivity. In some cases, such as molecular imaging for example, we know that the linear observations 
are formed as the Fourier transform of a positive measure. This introduces additional restrictions on the 
structure of these observations, which can be written as convex constraints on the phase vector. We detail 
two different ways of accounting for this positivity assumption. 

3.4.1. Direct nonnegativity constraints. In the case where the signal is nonnegative, [Waldspurger et al., 
2012] show that problem (2) can modified to specifically account for the fact that the signal is real, by 
writing it 

\\Ax - d\3i^{h)u\\l, 



mm 

^iGC^,|^i^|=l, 



using the operator T(-) defined as 



nz)^ 



Re(Z) - Im(Z) 
Im(Z) Re(Z) 

we can rewrite the phase problem on real valued signal as 



(7) 



mimnuze 



r{A) 



diag 



Re(ix) 



subject to i/ G C 



The optimal solution of the inner minimization problem in x is given by x = A\B2V, where 



A, 



Re(A) 
Im(A) 



B2 = diag 



and V 



Re{u) 
lm{u) 



hence the problem is finally rewritten 

minimize ||(742A^S2 - ^2)^112 

subject to vf + v^^- = 1, i = 1, . . . , n, 

in the variable v G M?^. This can be relaxed as above by the following problem 

minimize Tr{VM2) 

subject to Vi + V;\.^+. = 1, i = 1, . . . , n, 

which is a semidefinite program in the variable V G 82^, where 

M2 = (^24^2 - ^2)^(^24^2 - B2) = S|^(I - ^24)^2. 

Because x = A2B2V for real instances, we can add a nonnegativity constraint to this relaxation, using 

xx^ = {aIB2)uu^ {AlB2f 
and the relaxation becomes 

minimize Tr{VM2) 

subject to (4^2)1^(4^2)^ > 0, 

1/^0, 
which is a semidefinite program in y G S2n- 

3.4.2. Bochner's theorem. Another way to include nonnegativity constraints on the signal, which preserves 
some of the problem structure is to use Bochner's theorem. Recall that a function / : M^ ^^ C is positive 
semidefinite if and only if the matrix B with coefficients Bij = f{xi — Xj) is Hermitian positive semidefinite 
for any sequence Xi G W. Bochner's theorem then characterizes Fourier transforms of positive measures. 

Theorem 3.1. A function f : W \-^ C is positive semidefinite if and only if it is the Fourier transform of a 
(finite) nonnegative Borel measure. 

Proof. See [Berg et al., 1984] for example. ■ 

For simplicity, we first illustrate this in dimension one. Suppose that we observe the magnitude of the 
Fourier transform of a discrete nonnegative signal x G M^ so that 

\J^x\ — b 

with 6 G M^. Our objective now is to reconstruct a phase vector u e C^ such that |t^| = 1 and 

Tx — diag(6)ix. 

If we define the Toeplitz matrix 

Bijiy) = ^|i-j|+i, i,j = l,...,j), 
so that 

( yi yl •" yl\ 

y2 yi yl 



B{y) = 



2/2 yi 2/2 



y2 yi y2 
\ yn ... y2 yi J 

then when Tx = diag(6)i/, Bochner's theorem states that B{di3.g{b)u) ^ iff x > 0. The contraint 
B{di3.g{b)u) ^ is a linear matrix inequality in u, hence is convex. 
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Suppose that we observe multiple illuminations and that the k masks /i, . . . , //e G R^^^ are also nonneg- 
ative (e.g. random coherent illuminations), we have 

/ /io^\ 

V hoT J 

and the phase retrieval problem (5) for positive signals x is now written 

minimize u'Mu 

subject to Bj(diag(6)) ^0, j = 1, . . . , A; 
l^x^l = 1, z = 1, . . . n, 

where Bj (y) is the matrix B{y^^^) where y^^^ G C^ is the j^^ subvector of y (one for each of the k masks). 
We can then adapt the (PhaseCut) relaxation to incorporate the positivity requirement. In the one di- 
mensional case, using again the classical lifting argument in [Shor, 1987; Lovasz and Schrijver, 1991], it 
becomes 

min. Ty{UM) 

subject to diag(t/) = 1, i^i = 1, 

S^(diag(6)) ^0, j = 1, . . . , A; (PhaseCut+) 

U u 



u 1 

in the variables t/ G S^ and ^x G C^. The phase in u is fixed up to an arbitrary global shift, and the additional 
Til = allows us to exclude degenerate solutions with i/ = 0. 

Similar results apply in multiple dimensions, except that in this case the points xi in Bochner's theorem 
are located on a periodic grid. The matrix in this case is also significantly larger. 

4. Numerical Experiments 

We study molecular imaging problems based on electronic densities obtained from the Protein Data Bank 
[Berman et al., 2002]. From the 3D, we then obtain a 2D image by integrating the third dimension. After 
normalizing the image, we simulate diffraction observations through several masks. Here, masks consist of 
randomly generated binary filters. We simulate noisy observations the magnitude of the Fourier transform 
of the componentwise product of the image and the filter. As in the SPSIM package Maia [2013], we add 
random Poisson noise to the simulated observations, thus partly incorporating sensor noise. 

4.1. Parameters. We compare the results obtained by the Fienup algorithm and by the PhaseCut algorithm 
for the SDP relaxation while varying the number of masks and the noise level. For the PhaseCut relaxation, 
in order to deal with the large size of the lifted matrix, we use the low rank approximation described in §3.2 
and the sparsity of the observations vector described in §3.3. We then retrieve the phase vector as the first 
eigenvector in the final low rank approximation and refine it with the greedy algorithm Greedy or Fienup. 

More specifically, in the experiments that follow, the image was of size 128 x 128, we used a rank of two 
for the low rank approximation, kept the largest 1000 observations, did 5000 iterations of algorithm Fienup, 
and 20 cycles of algorithm PhaseCut (one cycle corresponds to optimizing once over all the columns of the 
lifted matrix). The seed of MATLAB random generator was set to 46. We compared the results of the phase 
recovery using 1 to 4 masks, and three different levels of Poisson noise (no noise, "small" noise, "large" 
noise). In all settings, all points of the electronic density were illuminated at least once by the random masks 
(the first mask lets all the signal go through). The noisy intensity measurements were obtained using the 
following formula, 

' max <; 0, a • Poisson 



where a is the input level of noise, and Poisson(A) is a random Poisson sample of mean A. Experiments 
were performed on a regular Linux desktop using Matlab for the greedy algorithms and a C implementation 
of the block coordinate algorithm for PhaseCut. CPU times are in seconds. 

4.2. Results. As shown in the table and figures above, in this setting, in most cases both algorithm Fienup 
and PhaseCut seem to converge to the (global) optimal solution, though Fienup is much faster. In some cases 
however, such as the experiment with two filters and no noise in Figure 2, initializing algorithm Fienup with 
the solution from PhaseCut significantly outperforms the solution obtained of algorithm Fienup which ap- 
pears to be stuck in a local minimum. The corresponding MSE values are listed in Table 2. In Figure 4 we 
plot the histogram of MSE for the noiseless case with only two illuminations, using either algorithm Fienup, 
or PhaseCut followed by greedy refinements, over many random illumination configurations. In many sam- 
ples, algorithm Fienup get stuck in a local optimum, while the SDP always converges to a global optimum. 
We also performa several simple experiments to illustrate the impact of positivity constraints on recon- 
struction performance. We generate an artificial signal x of dimension p with xi = i/p fori = 1, . . . ,p. 
We then test recovery performance using a varying number of random binary masks, for algorithms Fienup, 
the SDP relaxation PhaseCut with and without the positivity constraints described in §3.4 and the PhaseLift 
relaxation detailed in [Candes et al., 2011a], with and without nonnegativity constraints on the coefficients 
of the lifted variable X. The results are detailed in Table 1. 



Ilium. 


SNR 


Fienup 


SDP+ 


PhLift+ 


SDP 


PhLift 


1.00 


oo 


1.56 


1.21 


0.91 


1.31 


0.91 


2.00 


oo 


0.85 


0.56 


0.46 


1.75 


0.46 


3.00 


oo 


0.39 


0.02 


0.02 


1.52 


0.02 


4.00 


oo 


0.05 


0.03 


0.04 


0.75 


0.04 


1.00 


4.00 


1.67 


1.28 


0.91 


1.39 


0.91 


2.00 


4.00 


0.78 


0.76 


0.65 


2.21 


0.66 


3.00 


4.00 


0.40 


0.32 


0.27 


1.95 


0.31 


4.00 


4.00 


0.17 


0.19 


0.12 


1.39 


0.13 



Table 1. MSE over 10 samples in a random ID example for various values of the SNR, 
using a varying number of illuminations with random binary masks, for algorithms Fienup, 
the SDP relaxation PhaseCut with and without the positivity constraints (SDP+ & SDP) and 
the PhaseLift relaxation detailed in [Candes et al., 2011a], with and without nonnegativity 
constraints on the coefficients of the lifted variable X (PhLift+ & PhLift). 



5. User guide 

We provide here the instructions to artificially recover the image of a molecule from the Protein Data Bank 
using PhaseCut. This example is entirely reproduced with comments in the script testPhaseRecovery.m. 

5.1. Installation. Our toolbox works on all recent versions of MATLAB on Mac OS X, and on MATLAB 
versions anterior to 2008 on Linux. Installation only requires to put the toolbox folder and subdirectories on 
the Matlab path. Use for instance the command: 

>> addpath (genpath ( ' MYPATH/phaseRecoverylmagingToolbox' ) ) ; 

where MYPATH is the directory where you have copied the toolbox. 
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1 ill, a = 



2 ill., a = 



3 ill, a = 



4 ill, a^O 







1 ill., Q^ 




2 ill., Q^ 






1 ill., a 




2 ill., Q^ 




3 ill., a 




4 ill., Q^ 




Figure 2. Solution of the semidefinite relaxation algorithm PhaseCut followed by greedy 
refinements, for various values of the number of filters and noise level a. 



Nb masks 


a 


SDP obj 


SDP Refined obj 


Fienup obj 


SDP time 


Fienup time 


1 





0.410 


0.021 


0.046 


104 


78 


2 





2.157 


0.069 


0.877 


179 


156 


3 





4.036 


0.000 


0.000 


256 


241 


4 





7.075 


0.000 


0.000 


343 


343 


1 


10"^ 


0.806 


0.691 


0.704 


102 


67 


2 


10"^ 


3.300 


2.986 


2.992 


182 


144 


3 


10"^ 


5.691 


5.276 


5.277 


263 


229 


4 


10"^ 


8.074 


7.259 


7.259 


351 


308 


1 


10-^ 


3.299 


3.187 


3.193 


104 


71 


2 


10-^ 


7.492 


8.622 


8.646 


182 


145 


3 


10-^ 


10.419 


14.277 


14.288 


263 


222 


4 


10-^ 


14.139 


21.472 


21.444 


349 


305 



Table 2. Performance comparison between 
values of the number of filters and noise level 



algorithms Fienup and PhaseCut for various 
a. CPU times are in seconds. 



5.2. Generate the diffraction pattern of a molecule. Suppose we work with the caffeine molecule, on an 
image of size 128 x 128. We first load the molecule from PDB using the commands 

>> nameMol=' caf f eine .pdb' ; 
» N = 12 8 ; 
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1 ill., a = 



2 ill., a = 



3 ill., a^O 



4 ill., a^O 




• ^ 

- y* 





1 ill., a = 10" 



2 ill., a = 10" 



3 ill., a = 10" 



4 ill., a = 10" 







1 ill., a 




2 ill., a = 10" 



tesn 



3 ill., a 




4 ill., Q^ 




Figure 3. Solution of the greedy algorithm Fienup, for various values of the number of 
filters and noise level a. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

MSE 



Figure 4. Histogram of MSE for the noiseless case with only two illuminations, using 
either algorithm Fienup (red), or PhaseCut (blue) followed by greedy refinements, over 
many random starts. 



>> mol = pdbreadatom (nameMol) ; 

then calculate the electronic density rho of the molecule. Nz is the size of the density approximation (set 
here to 80), 



>> Nz 



80 
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>> [rho,Xr, Yr, Zr] = density (mol, [N, N, Nz ],' nbPts' ) ; 

Now, we set the parameters of the masks. The number of masks (also called filters or illuminations) is set 
to 2. Moreover we set the filter size to 64. The filter size is the square root of the number of pixels in the 
binary filter, Lcihc filter is of size filter size x filter size. The filter size must divide N (the square root 
of the number of pixels in the image). 

>> filter_size = 16 ; 
>> nb_f ilters=2 ; 

Since the filters are generated randomly, we set the seed of the uniform random generator to 1 in order to 
get reproducible experiments. Note that the quality of the phase retrieval may depend on the shape of the 
generated filters, especially when using only 2 or 3 filters. 

>> rand (' seed' , 1) ; 

Now we can generate an image, 2 filters and their corresponding diffraction patterns. We set the level of 
noise on the observations to zero here (i.e.no noise)., while a is the level of Poisson noise, and /3 is the level 
of Gaussian noise. 

>> alpha=0; 
>> beta=0; 

We set the oversampling parameter for the Fourier transform to 2. 

» M = 2 

The total number of observations, Le.ihc size of the vector b is 

>> nbObs=N^N^M^M^nb_f liters; 

Suppose that we want to use only the first largest thousand of observations in PhaseCut. We set: 

>> nbObsKept=1000; 

Note that the number of observations that is sufficient to get close to the optimal solution depends on the size 
of the data N and the sparsity of the vector b. From a more practical point of view, the larger nbObsKept, 
the more time intensive the optimization. Therefore, for a quick test we recommend setting nbObsKept to 
a few thousands, then increasing it if the results are not satisfying. 

Finally we call the function genData which is going to generate both the image x we want to recover, 
the filters, and the observations b. bs corresponds to the thousand largest observations, xs is the image 
recovered with the true phase but using only bs. idxijg is the logical indicator vector of bs (bs = b{idxi,s))' 
We put displayFig to 1 in order to display the filters, the images of the molecule x and xs, as well as the 
diffraction patterns (with and without noise). 

>> dlsplayFlg=l; 

>> [x, b, filters, bs, xs, ldx_bs] = genData (rho, Xr, Yr, Zr, .. . 

nb_fllters, fllter_slze, N, alpha, beta, Nz,M, nbObsKept , displayFig) ; 

5.3. Phase Retrieval using Fienup and/or PhaseCut. Using the data generated in the previous section, 
we retrieve the phase of the observations vector b. Suppose we want to use the SDP relaxation with greedy 
refinement, we set 

>> algorithm^' SDPRef Ined' ; 

The other choices for algorithm are Fienup, GreedylnPhase and SDP (no greedy refinement). Note that the 
greedy in phase algorithm is currently slower to run than Fienup (for equivalent results), so we rather use 
Fienup. 

We set the initial (full) phase vector u to the vector of ones, and the number of iterations for Fienup 
algorithm to 5000. The number of iterations for Fienup algorithm must be large enough so that the objective 
function converges to a stationary point. Usually 5000 iterations seem to be enough. We also set the number 
of iterations for the greedy in phase algorithm, although it will not affect our calculations. If we had chosen 
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the greedy in phase algorithm, one iteration would correspond to optimizing once each variable of the phase 
vector. 

>> uInit=ones (nbObs, 1) ; 

>> nbIterFienup=5000; 

>> nbIterGreedyInPhase=l; 

We also need to choose which algorithm we want to use in order to solve the SDP relaxation. We recommend 
to always use the block coordinate descent algorithm with a low rank approximation of the lifted matrix 
(BCDLR), since interior points methods (IP) and block coordinate descent without low rank approximation 
(BCD) become inefficient when the number of observations used is over a few thousands. 

» SDPsolver=' BCDLR' ; 

If we had wanted to use the IP solver or the BCD solver we would have set up respectively 

>> SDPsolver=' IP; 

or 

» SDPsolver='BCD' ; 

We can now set up the parameters for the BCDLR solver. 

>> nbCycles=20; 
» r=2; 

One cycle corresponds to optimizing over all the columns of the lifted matrix. In most cases, it seems 
that using nbCycles between 20 and 40 is enough to get close to the optimum, r is the rank for the 
low rank approximation of the lifted matrix. Similarly it seems that r between 2 and 4 gives reasonable 
results. Note that you can check that the low rank approximation is valid by looking at the maximum ratio 
between the last and the first eigen values throughout all iterations of the BCDLR algorithm. This ratio is 
outputted as eigRatio when calling the function retrievePhase (see below). We finally call the function 
retrievePhase in order to solve the SDP relaxation with greedy refinement. 

>> [retrievedPhase, objValues, finalOb j , eigRatio] ^retrievePhase (b, bs, .. . 
idx_bs, filters, M, algorithm, nblterFienup, nblterGreedylnPhase, . . . 
SDP solver, nbCycles, r, ulnit ) ; 

The function retrievePhase outputs the vector of retrieved phase u as retrievedPhase and the values of the 
objective function at each iteration/cycle of the algorithm in objValues. If using the SDP relaxation, the 
vector retrievedPhase is the first eigenvector of the final lifted matrix in PhaseCut. The "tic" "toe" enables 
to display the CPU time used to call the function retrievePhase. Note that the objective value in Fienup and 
in the SDP relaxation do not correspond exactly since the lifted matrix may be of rank bigger than 1 during 
the iterations of the BCDLR. Therefore we also output finalObj, which is the objective value of the phase 
vector extracted from the lifted matrix (i.e. the vector retrievedPhase). The image can now be retrieved 
using the command 
>> xRetreived=pseudo_inverse_A (retrievedPhase . ^b, filters, M) ; 

Finally you can visualize the results using the following standard Matlab commands, plotting the objective 
values 

>> figure (3) 

» subplot (2, 1, 1) ; 

>> title (algorithm) 

>> plot (loglO (abs (objValues) )) ; axis tight 

and displaying images 

» subplot (2, 3, 4) 

>> imagesc (abs (x) ); axis off; 

>> subplot (2, 3, 5) 
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>> imagesc (abs (xs) ) ; axis off; 

>> subplot (2, 3, 6) 

>> imagesc (abs (xRetreived) ) ; axis of f ; 

5.4. Reproduce the experiments of the paper. All the numerical experiments of the paper can be repro- 
duced using the Matlab scripts included in the toolbox. Use the script testSeeds.m to generate the 
histogram in Figure 4. Use the script testNoiseNblllum . m to generate Figure 2, Figure 7 and Table 2. 
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